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ABSTRACT 

L^:fJl P iT a ' C ° ntr0lpr °?i em arising in non-planar orbital transfer employing 
enprnfS IS addressesd - The mission involves the transfer from high 

QtraSnw 0 ? * t0 ^ energy orblt ( LE0 ) with orbital plane change. The basic 
strategy here is to employ a combination of propulsive maneuvers in space and 

J^T?Pn n t the f tmos P here - The basic sequence of events for the 
H u° ° LE T? transfer consists of three phases. In the first phase, the 
orb'tal transfer begins with a deorbit impulse at HEO which injects the vehicle into an 

!oh P i^ tranSfe t r ° r ^ With perigee inside the atmosphere. In the second phase the 

riPk^H prh? P in J a y c ° ntrolled by lift and bank angle modulations to perform the 
deisred orbital plane change and to satisfy heating constraints. Because of the energy 

°p S H l U SJ ha u m ’, an ,rap , u ' sa is squired to initiate the third phase to boost the 
vehicle back to the desired LEO orbital altitude. The third impulse is then used to 

wNrh a | l |cti he orblt at LE 9' The problem is solved by a direct optimization technique 
onH o U P'ecewise polynomial representation for the state and control variables 

nn«nS 0Ca l ,0r I t0 tf sfy - the differential equations. This technique converts the 
nnmTrv ‘T i 0011 ? . P * r0b em mt0 a nonl,nea r programming problem which is solved 
™ C i' y ' ( Solutions were obtained for cases with and without heat constraints and 
tor cases of different orbital inclination changes. The method appears to be more 

hpnHip U and i r °^ USt than , other optimization methods. In addition, the method can 
handle complex dynamical constraints. 
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NOMENCLATURE 


A : Sps/2 

Cd : drag coefficient 

Cdo : zero-lift drag coefficient 

Cl : lift coefficient 

Clr : lift coefficient for maximum lift-to-drag rati 

D : drag force 

g : gravitational acceleration 

g s : gravitational acceleration at surface level 

H : altitude 

J : performance index 

K : induced drag factor 

L : lift force 

m : vehicle mass 

R : distance from Earth center to vehicle center of gravity 

Ra : radius of the atmospheric boundary 

Rc : radius of the low Earth orbit (LEO) 

R d : radius of the high Earth orbit (HEO) 

Re radius of Earth 

S aerodynamic reference area 

t : time 

V : velocity 

T : Thrust 

P : inverse atmospheric scale height 

y : flight path angle 

\y : heading angle 

a : bank angle 

0 : down range angle or longitude 

<]> : cross range angle or latitude 

p : gravitational constant of Earth 

p : density 

AV : characteristic velocity 

subscripts 

c : subscript for circularization or reorbit 

d : subscript for deorbit 

s : subscript for surface level 


1. INTRODUCTION 

In order to have a viable and affordable space program, advanced technology must be 
exploited and new design concepts must be developed to reduce the size and cost of 
transportation elements for supporting new mission requirements. One of the new 
concepts that has evolved in recent years to advance the cost effectiveness of space 
transportation systems is the aerodynamically assisted orbit transfer. Such an orbital 
transfer vehicle is designed with an aerodynamic configuration which can utilize the 
planetary atmosphere for the purpose of energy management. Numerous studies 
have demonstrated that the use of the aerobraking can significantly reduce the 
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propulsive velocity requirements for certain class of orbit transfers. Excellent review 
papers were given by Warberg (Reference 1) and Mease (Reference 10). 

In this paper, the fuel optimal control problem arising in a typical orbit transfer with 
plane change from HEO to LEO as discussed in most recent publications is 
addressed. In this case, as discussed in Reference 2, the aeroassisted orbit transfer 
vehicle (AOTV) maneuver involves three propulsive burns or impulses as sketched in 
Fig.1 . In the first phase , the transfer begins with a deorbit impulse at HEO which 
injects the vehicle into an elliptic transfer orbit with the perigee inside the atmosphere. 
In the second phase, the vehicle is inside the atmosphere and is optimally controlled 
by the lift and bank angle modulations to perform the desired orbital plane change and 
to satisfy the heating rate and other physical constraints. Because of the the energy 
loss during the atmospheric maneuvers, an impulse is required to initiate the third 
phase to boost the vehicle back to the final orbital altitude. Finally, the third impulse is 
applied to circularize the orbit at LEO. In summary, there are three propulsive burns 
and an aeroassist plane change inside the atmosphere. Simulation results similar to 
those obtained in the draft paper of Reference 2 have been obtained here by using the 
Hermite polynomial and collocation technique to convert the optimal control problem 
into a nonlinear programming ( NP ) problem which is solved numerically using the 
optimization code, NZSOL ( cf. Reference 12 ) provided by Gill, which is an improved 
version of NPSOL ( cf. Reference 6 ), developed at Stanford. This solution method is 
different from the indirect method such as those discussed in Reference 2,4,7 and 8. 
The above simulation results have been extended to cases with heating constraints 
and cases for different orbit inclination changes. The details are presented and 
discussed in this paper. It is important that in the future these simulations be extended 
to include all other realistic flight constraints and to establish baseline optimum 
trajectory characteristics for GEO to Space station or shuttle, lunar and Mars missions. 


2. DIRECT TRAJECTORY OPTIMIZATION WITH COLLOCATION AND 
HERMITE POLYNOMIALS 

In the direct collocation with nonlinear programming approach, the trajectory is 
approximated by piecewise polynomials, which represent the state and control 
variables at a number of discrete time points, i.e., nodes. For a given state variable, 
the state trajectory over a given "segment" between two nodes is taken to be the 
unique Hermite cubic which goes through the end points of the segments with the 
appropriate derivatives that are dictated by the differential equations of motion at the 
endpoints. This is the "Hermite cubic" since it is determined by the states and their 
derivatives. A collocation is taken at the center of the segment where the derivative 
given by the Hermite cubic is compared to the derivative obtained from the evaluation 
of the equations of motion. The difference is termed the "defect" and is a measure of 
how well the equations of motion are satisfied over the segments. If all the defects are 
zero, then the differential equations are satisfied at the center collocation points as 
well as at the endpoints. Figure 2 shows the typical defects between node 1 and node 
2 . 
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(2-1 a) 


Let the system of equations of motion be given as 


X’ = f(X,U) 

where X is the state vector and U is the control vector and (') denotes the differentiation 
with respect to the time. Let the time over a given segment be T. For the problem 
discussed here, one can show that 


X = ( x, y, z, x’, y\ z\ m ) 
U-(C Ll o) 


(2- lb) 


Then the Hermite interpolated x-component of the state vector X at the center point is 
x c = (1/2) (X! + x r ) + (T/8) [f(X 1t Ui) - f(X r ,U r )] (2-2) 


where xi and x r are respectively the x-component of the state vector X at the left and 
the right nodes. The derivative of the interpolating Hermite cubic at the center point is 

x c ' = -3/(2T) (x, • x r ) - (1/4) [f(X,,Ui)+ f(X r ,U r )l (2-3) 


The defect vector is then calculated as 


d = f(X c ,U C ) - x c ’ 


(2-4) 


If xi, ui, x r , and u r are chosen such that the elements of the defect vector, d, are 
sufficiently small, the "Hermite polynomials" become an accurate approximation to the 
solution of the differential equations of motion (by implicit integration). With the above 
approach, the differential equations are converted into nonlinear ajgebraic equations 
and the optimal control problem can then be solved using the nonlinear programming 
techques. 

3. APPLICATION TO OPTIMAL AEROASSISTED ORBITAL TRANSFER 
WITH PLANE CHANGE 

The aeroassited orbital transfer can be analyzed in three phases , i.e., deorbit, 
aeroassist (or atmospheric flight), boost and reorbit (or circularization). In each of the 
phases, a particular set of equations of motion apply. 

3.1 Deorbit 

Initially, the spacecraft is moving with a circular velocity Vd = V^/Rd ' n a circular orbit 
of radius Rd, well outside the Earth’s atmosphere. Deorbit is accomplished at point D 
by means of an impulse AVd, to transfer the vehicle from a circular orbit to an elliptic 
orbit with perigee low enough for the trajectory to intersect the dense part of the 
atmosphere. Since the elliptic velocity at D is less than the circular velocity at D, the 
impulse AVd is executed so as to oppose the circular Velocity Vp . In other words, at 
point D, the velocity required to put the vehicle into elliptic orbit is less than the velocity 
required to maintain it in circular orbit. The deorbit impulse AVd causes the vehicle to 
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enter the atmosphere at radius R a with a velocity V e and flight path angle y e . It is 
known that the optimal energy loss maneuver from the circular orbit is simply the 
Hohmann transfer and the impulse is parallel and opposite to the instantaneous 
velocity vector. 

After applying the deorbit impulse and before entering the atmosphere at R a , the 
deorbit trajectory is a coasting arc and known integrals of the equations of motion can 
be used to relate the state vectors at R a ,the entry into atmosphere to the state vectors 
right after the deorbit impulse at Rd. Using the principle of conservation of energy and 
angular momentum at the deorbit point D and the atmospheric entry point E, we get 

V 2 /2-ji/R a =(V d -AV d ) 2 /2-n/R d ,3.,, 

R a V e cos(-y e ) = R d (V d - AV d ) 


from which we can solve for AVd to get 

AV d = 1 R d - ^2p(1 / R a - 1 / R d )/[(R d / R a ) 2 / cos 2 y e - 1] 

It is easily seen that the minimum deorbit impulse AVdm obtained for y e = 0, 
corresponds to an ideal transfer with the space vehicle grazing the atmospheric 
boundary. To ensure proper atmospheric entry, the deorbit impulse AV d must be 
higher than the following minimum deorbit impulse AVdm 

AV dm = Vn / Rd " ^(1 ' R a - ' 1 ■ 7 \ Rd) / [(R d / R a ) 2 - 1] 

(3-4) 

Physically, the second term of the above equation corresponds to the apogee velocity 
of an ellpitic transfer orbit with perigee radius R a and apogee radius Rd. This elliptic 
transfer orbit is tangent to the atmosphere boundary at perigee. It will be shown later 
that the nonlinear constraint equations ( 3-15 ) at the atmospheric entry point can also 
be derived from equations ( 3-1 and 2 ). 

3.2 Aeroassist 

During the atmospheric flight, the vehicle is optimally controlled by the lift and bank 
angle modulations to achieve the necessary velocity reduction (due to the atmospheric 
drag) and the plane change. In the present formulation, only the aeroassisted 
atmospheric flight need be solved by using the collocation and nonlinear 
programming techniques discussed earlier in this paper. The solutions in the other 
phases are provided by the known integral relations of the equations of motion 
because these arcs are coasting arcs. 

Consider a vehicle with the point mass m, moving about a rotating spherical planet. 

The atmosphere surrounding the planet is assumed to be at rest, and the central 
gravitational field obeys the usual inverse square law. The equations of motion for the 
vehicle are given by (Figure 1), 
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r = Vsiny 

(3-5a) 

^ V cosy cosy 
rcos<]> 

(3-5b) 

1 

• Vcosysiny 
9 = 

r 

(3-5c) 

v (TiTcose-D) psiny +C0 2 rcos< j) ( s i n y C0S <j> — eosy siny si n< » 
m r 2 

(3-5d) 

(r|T sine + L)cosa p cosy + Vcosy rns^ 

r mV V r 2 r 

2 

+ (D rc0S< l > ^ COS y cos9 + siny S jn\(/ sin<J>) 

i 

(3-5e) 

(tlT sine+L)sino V cosy cosy tan$ + 2(j)(tanY siny cos0 _ sin(rt 
y mV cosy r 

to 2 r cos y sin9 cos 9 
V cosy 

(3-5f) 

m = -f(r,V,Ti) 

(3-5g) 

where for a given vehicle, the drag D and the lift L are 

D = ^-pV 2 C D 
2m 

(3-5h) 

L = A p v 2 C L 
2m 

(3-51) 

and the drag and lift coefficients obey the drag-polar relation 
Cd = Cdq + KCl 

(3-Sj) 

Also, for an exponential atmosphere, one has 
p = p s exp(-Hp) and H = R-R E 

(3-5k) 


Simulation results obtained here were using the U.S. standard Atmosphere 1976. 
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Perth® problem considered here, one assumes that, inside the atmosphere, the 
vehicle is optimally controlled by the aerodynamic forces only. Thus it is assumed that 
the thrust T is absent and the point mass is constant in this region. Furthermore no 
earth rotation was assumed. The later is equivalent to consider the motion with 
respect to an earth fixed inertial coordinate system ( ECI ). The plane change or the 
orbit inclination, /, is related to the cross range <|> and the heading angle \jr as 


cosi = cos <|> cosy t e <t<t f 


The orbit inclination changes throughout the atmospheric flight and must end up with 
the required value at exit. For small values of cross range angle O, /' is given by the 
heading angle \j / itself. 


3.3 Boost and Reorbit 

Ourmg the atmospheric flight, the vehicle undergoes the plane change using the lift 
and bank angle modulation . Because of th© loss of energy during the atmospheric 
maneuver, a second impulse is required at the exit from the atmosphere to boost the 
vehicle back to the final orbital altitude at LEO. 


The vehicle exits the atmosphere at point F, with a velocity Vf and the flight path angle 
Yf. The additional impulse AV b , required at the exit point F for boosting the vehicle 
into an elliptic transfer orbit with apogee radius R c> and the reorbit (or circularization) 
impulse A V c , required to insert the vehicle into a circular orbit, are obtained by using 
the principle of conservation of energy and angular momentum at the exit point Fand 
the reorbit or circularization point C. Thus, we have 


(Vf + AV^) 2 / 2 - ji / R a = (V c - AV C ) 2 / 2 - |i / R c 
(V, + AV b )R a cosy, = R C (V C - AV C ) 

Solving for A V b and AV C from the above equations (3-7) and (3-8) yields 


AV b = J 2p(l / R a - 1 / R c ) / 1 - (R a / R c ) 2 cos 2 Yf 1 - V, 

Y L 1 (3-9) 

AV C = Jn / R c - J 2 p(1 / R a - 1 / R c ) / [(R c / R a ) 2 / cos 2 Y f - ll 

v 1 J (3-10) 

It is interesting to note that the second term of equation (3-10) is maximum for y =0 
and therefore the reorbit impulse AV C ,is minimum for Yf =0. It will be shown later that 
boundaiy conditions and nonlinear constraint equations at the exit point F, can be 
denved in terms of the final orbit characteristics and the final state vectors at the exit as 
shown in (3-1 6,17, & 18). 
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3.4 Performance Index 

It is known that the change in speed, AV, also called the characteristic velocity, is a 
convenient parameter to measure the fuel consumption. For minimum-fuel maneuver, 
the objective is then to minimize the total characteristic velocity. A convenient 
performance index is the sum of the characteristic velocities for deorbit, boost, and 

reorbit, as 

J = AV d + AV b +AV c (3-11) 

Where, AV, AVt>, and AVq are the deorbit, boost, and reorbit characteristic velocities 
respectively, and are related as 


AV d = VP/Rd “ ( R a / R d) v e COS(-Ye) 


(3-12) 


AV C = >/R c -(R a / R C )(V| + AV b )cosy f (3-13) 

Alternatively, AVd, AVt* and AV C are also given by (3-3, 9, and 10) respectively. Note 
that for a given final circular orbit, the impulses AVb and AVq are completely 
determined by the state variables Vf and Yf at the exit of the atmospheric portion of the 
trajectory. The velocity V e and the flight path angle Ye at the atmospheric entry point 
are dependent only on the magnitude of the deorbit impulse AVd . It follows that the 
optimal control problem needs to consider only the trajectory segment within the 
atmosphere subject to the nonlinear constraints and boundary conditions at the 
atmospheric entry and exit points. In addition, other path constraints such as the peak 
heating rate have to be satisfied. 

3.5 Boundary conditions and constraints 

The boundary conditions and constraints for the optimal control problem can be 
summarized as follows: 

• At the entry into atmosphere, the following initial constraints. must be satisfied. 

R = R a ; Ye^° • <t>e =°- Ve = 0 (3-14) 


Ve 2 

1- 

— 1 cos 2 (Ye) 

-p 

f— 1 

2 


l R dJ 


< R a R d, 


The first initial constraint is required to ensure the vehicle enters the atmosphere. In 
the present formulation, the initial velocity and the flight path angle are unknown and 
to be determined by the optimization processes subject to the second constraint. 

• At the exit from atmosphere, the following constraints must be satisfied. 


R = R a : Yt -0 


(3-16) 
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(3-17) 





COS if - COS <t>f COS \l/f = 0 

Vf (3-18) 

Equation (3-16) is required to ensure the vehicle exit the atmophere. The second 
costraint can be used to compute AVb, and if AVb , is assumed to be zero as in the case 
of aerobraking without orbit plane change, the above constraint must be imposed to 
deternm ne the correct Vf and The third constraint is required to perform the desired 
orbital plane change. 

• In addition, there are other path constraints ,i.e., constraints must be satisfied 
along the trajectory such as 

a) Stagnation Point Heating Rate Constraints 

- b) Altitude Constraints 

- c) Bounds on the Control Variables 

- d) Others 

4 ‘ iH2Ji? TURE AND SOLUTION OF THE NONLINEAR PROGRAMMING 
PROBLEM 

The direct collocation and Hermite polynomial procedures described above convert 
optimal control problems into corresponding nonlinear programming problems. 
Ordinary differential equations are converted into corresponding nonlinear algebraic 
equations (or nonlinear “defects” constraint equations). These problems can then be 
solved using nonlinear programming codes. 

The variables for the nonlinear programming problem are the collected state vectors 
and control vectors at the nodes and the time duration of phases. These quantities are 
assembled into the NLP state vectors 

* T= [x7>Ui\ Xn,Un,ti,t2 *kl 

J (4-1) 

where n is the number of nodes and k is the number of phases on the trajectory. The 
defects and other physical and mathematical constraints are collected into the NLP 
constraint vector C 

cT = [ d l < d 2 > dj, w[, wj.wj. w7] 

L 1 J (4-2) 

where dj is the defect vector and w is a vector of additional problem constraints. 

The nonlinear programming code used here is the NZSOL (Reference 12). The 
NZSOL is an improved version of the NPSOL (Reference 6) , developed by the 
Stanford Optimization Laboratory and designed to minimize a smooth nonlinear 
function subject to a set of constraints which may include simple bounds on the 
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variables, linear constraints, and smooth nonlinear constraints The problem is 
assumed to be stated in the following form: 


NP 


minimize F(x) 
x e R n 


subject to t < 


&)j 


s u, 


(4-3) 


where the objective function F (z) is a nonlinear function, Al is an mi_, x n constant 
matrix of general linear constraints, and c(x) is an m/v - vector of nonlinear constraint 
functions, the objective function F and the constraint functions are assumed to be 
smooth, i.e., at least twice-continuously differentiable. (The method of NPSOL will 
usually solve NP if there are only isolated discontinuities away from the solution). 

Note that upper and lower bounds are specified for all the variables and for all the 
constraints. This form allows full generality in specifying other types of constraints. In 
particular, the i-th constraint may be defined as an equality by setting £\ = uj. If certain 
bounds are not present, the associated elements of t or u can be set to special values 
that will be treated as - 00 or + . 


Here we briefly summarize the main features of the method of NZSOL and NPSOL as 
discussed in Reference 6 because Reference 12 is not available to general public. At 
a solution of NP, some of the constraints will be active, i.e., satisfied exactly. An active 
simple bound constraint implies that the correspondin g variabl e is fixed at it? boun d, 
and hence the variables are partitioned into fixed and free variables. Let C denote 
the m X n matrix of gradients of the active general linear and nonlinear constraints. 
The number of fixed variables will be denoted by npx, with hfr (hfr = n - hfx) the 
number of free variables. The subscripts “ FX ” and “FR” on a vector or matnx will 
denote the vector or matrix composed of the components corresponding to fixed or 
free variables. The details are discussed in Reference 1 1 . 


A point x is a first-order Kuhn-Tucker point for NP if the following conditions hold: 


(i) x is feasible; . . , , 

(ii) there exist vectors £ and X (the Lagrange multiplier vectors for the 

bound and general constraints) such that 


g = C T X + C, (4-4a) 

where g is the gradient of F evaluated at x, and = 0 if the j-th variable 
is free. 
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(m) The Lagrange multiplier corresponding to an inequality constraint 
active at its lower bound must be non-negative, and non-positive for 
an inequality constraint active at its upper bound. 


Let Z denote a matrix whose columns form a basis for the set of vectors orthogonal to 

the rows of C F r\ i.e., C F rZ=0. An equivalent statement of the condition in terms of Z 
is 


Z T 9fr = 0 


(4-4b) 


The vector ZTg FR j s termed the projected gradient of F at x. Certain additional 

conditions must be satisfied in order for a first-order Kuhn-Tucker point to be a solution 
of NP. 


4-1 The Quadratic Programming Subprobiem 

Similar to NPSOL, the basic structure of NZSOL involves major and minor iterations. 
The major iterations generate a sequence of iterates ( x i<) that converge to x\ a first- 

order Kuhn-Tucker point of NP . At a typical major iteration, the new iterate x is 
defined by 


where x is the current iterate, the non-negative scalar a is the step length, and p is the 
search direction. Also associated with each major iteration are estimates of the 
Lagrange multipliers and a prediction of the active set. 


The search direction p is the solution of a quadratic programming subproblem of the 
form 

minimize gVVnp 


subject to 


l< 


P 

A L p 

[AnP 




(4-5b) 


where g is the gradient of F at x, the matrix H is a positive-definite quasi-Newton 
approximation to the Hessian of the Lagrangian function and An is the Jacobian 
matrix of c evaluated at x. 


The estimated Lagrange multipliers at each major iteration are the Lagrange 
multipliers from the subproblem (and similarly for the predicted active set) and provide 
information about the the sensitivity of these NLP problems. 

Certain matrices associated with the QP subproblem are relevant in the major 
iterations. Let the subscripts "FX” and “FFt” refer to the predicted fixed and free 
variables, and let C denote the m x n matrix of gradients of the general linear and 
nonlinear constraints in the predicted active set. First, we have available the TQ 
factorization (Reference 1 1 ) of C F r : 
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Cfr Qfr = (0 T ), (4-6) 

where 7 is a nonsingular m x m reverse-triangular matrix (i.e., f,y = 0 if / +j < m), and 
the non-singular dfr x hfr matrix Qfr is the product of orthogonal transformations. 
Second, we have the upper-triangular Cholesky factor R of the transformed and re- 
ordered Hessian matrix 


R t R = H q = Q t HQ, (4-7) 

where H is the Hessian H with rows and columns permuted so that the free variables 
are first, and Q is the n x n matrix 


( Q, 


Q = 


! fr 


\ 


•fx. 


(4-8) 


with Ifx the identity matrix of order hfx- If the columns of Qfr are partitioned so that 


Qfr = ( z Y), (4-9) 

the riz ( n z = n FR~ m ) columns of Z form a basis for the null space of Cfr • The matrix 
Z is used to compute the projected gradient Z T gFR at the current iterate. 


As discussed in Reference 6 and 1 1 , a theoretical characteristic of SQP methods is 
that the predicted active set from the QP subproblem is identical to the correct active 
set in a neighborhood of x*. In NPSOL, this feature is exploited by using the QP active 
set from the previous iteration as a prediction of the active set for the next QP 
subproblem, which leads in practice to optimality of the subproblems in only one 
iteration as the solution is approached. Separate treatment of bound and linear 
constraints in NPSOL also saves computation in factorizing Cfr and Hq. 


4.2 The merit function . _ . ... . 170r ., . 

Detailed discussions of the merit function are given in Reference 14. In NZbUL and 
NPSOL, once the search direction p has been computed, the major iteration proceeds 
by determining a steplength a that produces a “sufficient decrease in the augmented 
Lagrangian merit function 


L(x,X,s) = F(x)-X Xj(Cj(x)-Si)+^S PifciM-Si) 2 , 

i * i (4-1 U) 

where x, X and s vary during the line search. The summation terms involve only ths 
nonlinear constraints . The vector X is an estimate of the Lagrange multipliers for the 

nonlinear constraints of NP. The non-negative slack variable (Sj) allow nonlinear 
inequality constraints to be treated without introducing discontinuities. The solution of 
the QP subproblem (4-5) provides a vector triple that serves as a direction of search for 
the three sets of variables. 


4.3 The quasi-Newton updated 
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goin 9 int0 the deta iled discussions, it is important to point out that both the 
NZSOL and NPSOL start by initializing the Hessian matrix H = Identity matrix. Thus at 
the beginning, the search direction is in the steepest decent direction. No initial 
curvature information is computed and the curvature information is accumulated 
through the BFGS quasi-Newton updates. The matrix H in (4-5) is a positive-definite 
quasi-Newton approximation to the Hessian of the Lagrangian function. At the end of 
each major iteration, a new Hessian approximation H is defined as a rank-two 
modification of H. In NPSOL the BFGS quasi-Newton update is used: 


H = H- 


-y HSS T H + 

s'Hs 



where s = x-x (the change in x ). 


(4-11) 


Rather than modifying H itself, the Cholesky factor of the transformed Hessian Hq (4- 
7) is updated, where Q is the matrix from (4-8) associated with the active set of the QP 
subproblem. The update (4-1 1) is equivalent to the following update to Hq : 


*j 

Hq = h q ~ rr H q h q s q s q Hq + ~ T ~ y Q yQi 

SqHqSq YqSq ^ 

where y Q = Q T y and s q = Q t s. This update mav be ex pressed as a rank-one 

tipcat? to B. and is used to incorporate new curvature information obtained in the 
move from x to x . 


4.4 NZSOL, NPSOL 4.02, and NPSOL 2.1 

For those who are interested in applying these NLP codes, there are two publised 
versions of NPSOL. The NPSOL 4.02 was developed after the NPSOL 2.1 and 
therefore more reliable and efficient algorithm were incorporated according to Gill ( 
Reference 12 ). However, in updating the Cholesky factor, the NPSOL 4.02 updates 
the whole or complete R while the NPSOL 2.1 updates only the part associated with 
Ihe Z-space or null space of R. For the problem formulated here .usually several 
hundred varibles are involved and the NPSOL 2.1 converges in less computing time. 
The NZSOL (Reference 12) incorporates not only latest efficient and reliable algorithm 
but also updates only the part of R associated with the null space of R only. In addition 
to improve the algorithm of NPSOL, it also adopts the best parts of both NPSOL 2 1 
and 4.02. 


Finally, it may be interesting to point out that the matrices in the present formulation 
using collocation and Hermite polynomial are large and fairly sparse. For 
computational efficiency, it is important to incorporate NLP codes such as MINOS 
(Reference 1 3) to take advantage of the special characteristic of the collocation 
formulation discussed here. 

5. NUMERICAL RESULTS AND DATA 

The data used in the numerical experiments presented here (c.f. Reference 2 and 9) 
are summarized as follows: 


Cqo = 0.1 ; K = 1.111 ; m/S = 300kg/m2 


(5-1) 
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and the drag polar is 


C D = C D o + K * C L 2 

(5-2) 

and other data are 



p a = 1.225 kg / m 3 ; \i = 3.986x1 0 14 m 3 / sec 2 
p = 1/6900 m" 1 ; R E = 6378km 

H a = 1 20 km; R d = 1 2996 km; R c = 6558 km ( 5-3 ) 

Using the above mentioned data, simulations were carried out. The optimal solution 
for the reference case ( shown in figures as Case 1 ) has the following entry and exit 

status. 


Entry status: H e = 1 20 km; V e = 9034.74 m/sec 

Ye = -4.36 degrees; <t> e = 0; y e = 0 ( 5-4 ) 

Exit status: H f = 1 20 km; V f = 7028.95 m/sec 

y f = 0.0 deg; = -6.69 deg 

Y f = 1 8.891 deg; total flight time = 478 sec ( 5-5 ) 


Characteristic velocities: 

Deorbit characteristic velocity 
Boost characteristic velocity 
Reorbit characteristic velocity 
Total characteristic velocity 


AVd = 1 031 .59 m/sec 

AVb = 821 .49 m/sec 

AV C = 1 7.98 m/sec 

AV = 1 871 .07 m/sec ( 5-6 ) 


Time histories of altitude, velocity, flight path angles, heading angles, dynamical 
pressure, atmospheric density, orbit inclination and heating rate are shown in Figure 
3-10. Figures 11-13 show lift coefficient, bank angle, and lift to drag ratio as a function 
time for several simulation runs for the reference case(.i.e., Case 1 ). These simulation 
runs show that at high altitude the control may be different for different simulation runs 
depending upon initial guesses. This is really not a surprise because at high altitude, 
the aerodynamic forces or the controls are ineffective. In fact, the problem has a weak 
optimium with respect to the controls at high altitude. 

Without going into the details, the characteristic velocities for the cases with orbital 
inclination changes for 15, 20 and 25 degrees are summarized in Table 2. 

The heating rate Q r t along the atmospheric trajectory, is computed for the stagnation 
point of a sphere of radius of one meter, according to the following relation 
( Reference 2 and 6 ) 

Qr-K r p0.5v3.08 ( 5 ' 7 } 
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i^ h thL e r!r^i® the r f tmos P heri c density in kg/km3, V is the velocity in km/sec and the K r 
is the proportionality constant equal to 0.000308. The peak heating rate for the 

reference case is about 239 W/cm2 Simulation results for cases with peak heat rate 

3 i n oTrfrnmnLk°n 3 n an ? 170W/cm2 were also obtained and shown in Figures 

t?. 0 w ! th the reference case without heat constraints as Case 2 and 3 

bv fs and q'n r L h r^ e + tW0 CaS ? S '! edu ^. es the P eaki ng heating rate of the reference case 
cLe*i?? d 3 .° P® rc f nt respectively. Simulation results presented here provide the 
sensitivity of trajectory and associated physical variables as the heating constraints 
die imposed. 

Similarly, the characteristic velocities for the cases with heating constraints are 
^, a .H Zed i" Table 2 - The u percent reductions are with respect to the peak heating 
thlrmlfo ? fe ; ence case , wlthout heat constraints. As shown here, one needs less 9 
thermal protection materials and more fuel consumption to fly the heat constrained 

S°[' eS and the [. ef ° re b V takin 9 ir| to account the weight of thermal protection 
materials one may find an optimal design to minimize the total vehicle weight. 

Another interesting observation from the data given in Table 1 and 2 is that the deorbit 

I/oKi t 6 f S a m ? st the s . ame for al1 the cases simulated here. The total characteristic 
velocity for a given optimal trajectory is almost completely determined by the boost and 

^circuiatfon. In fact, the boost velocity contributes the most to the variation of the 

it SI ^ tenStlC veloc,ty - Ph y s,cal| y- n is obvious as the vehicle makes a larger turn 

orbital SndT°Sit^^ needs more velocit y to boost it back to the final 

nt - d . A t *? ugh t . he tota ! chara cteristic velocity is insentive to the magnitude 
of deorbit impulse, the optimal trajectory is very sensitive to AV^ . 

6. CONCLUDING REMARKS 

An excellent survey of the subject was given in Reference 1. Walberg reviewed the 
problem of synergetic plane change for optimal orbital transfer. In a recent paper by 
Naidu (c.f. Reference 2), fuel optimal trajectories of aeroassisted orbital transfer with 
plane change were presented using the so-called multiple shooting method for the 
case without heat constraints and under the assumption that all the synergetic plane 
change was performed entirely in the atmosphere. A brief review of the progress 
made in this field was also given in Reference 2. In our paper, a similar problem for 
cases with and without stagnation point heating rate constraints was solved using the 
collocation and nonlinear programming technique. This method is especially suitable 
tor parametrical studies because of its relative insensitivity to initial guesses. Once a 
solution for a reference case is obtained, solutions for other cases such as different 
orpital inclination change can be obtained easily. 

Finally, the present problem can also be formulated under a more general assumption 
"f , n h ot al ‘ Ij® P |ane changes are entirely made in the atmosphere. It must be noted 
mat tne AOTV transfer can be made more efficient propulsively if the plane change is 
performed partly in the atmosphere and partly in space and the propulsive plane 
cnange in space is subdivided into components associated with various impulsive 
pomts. For the more general formulation discussed, the desired plane chanqe mav 
consist of more than one plane change, i.e. 

Total orbital plane change = /y + / 2 + '13 + i a 
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where ii h . and 13 are plane changes at the deorbit, reboost, and reorbit 
respectively and i a is the aeroassisted plane change and all these plane changes will 
be determined by the optimization processes discussed here. Preliminary simulation 
results were obtained and are to be published in the near future. 


It should be mentioned that the collocation and nonlinear programming technique 
discussed here was recently applied to another group of orbital transfer problem by 
Enright and Conway in Reference 3 and the relative insensitivity of this method to the 
initial guesses was also observed by them. Our basic simulation test bed is the OTIS 
codes ( Reference 5 ) with an improved and updated nonlinear programming code ( 
NZSOL ). All physical models used were documented in Reference 5. Of course, 
necessary modifications and corrections have to be incorporated to simulate the 
aerobraking problems discussed here. 


It may be worthwhile mentioning that the present problem was actually solved by 
guessing the initial state and control varibles at four selected points , i.e., the initial 
point, the final point and two other nodal points along the trajectory inside the 
atmosphere. The initial state and control variables at other nodes or grid points were 
simply obtained by linear interpolation. These initial guesses do not have to satisfy 
either the governing equations or the nonlinear constriants including the defects. Only 
roughly guesses are needed at these four points. Converged solutions were obtained 
with relative ease. However, it is important to point out proper scaling of the defects, 
constraints and variables are essential to get converged solutions. Although our 
results were compared with the draft paper of Reference 2, the solution presented by 
Naidu was not actually optimal because the final flight path angle yf = -0.621 7 degrees 
is negative. For simulations discussed here, converged solutions were obtained by 
using as little as 20 nodes. However, in some cases, converged solutions were 
obtained using 60 nodes, In the later case, the problem has more than 660 
independent variables and more than 400 nonlinear “defects" equations. For cases 
with heating rate constraints, the problem has more than 500 nonlinear constraint 
equations includinding the “defect" equations. As far as we know, this may be the first 
time converged soultions were obtained for so many independent variables and 
nonlinear constraint equations. This also illustrates how powerful the nonlinear 
programming code and the collocation and Hermite polynomial technique are. 


Finally, it is important to mention again that aeroassisted orbital transfer introduces a 
strong coupling between the vehicle design and the trajectory design as indicated by 
the simulation data. A trajectory that minimizes fuel mass, without attention to heating , 
may require the vehicle to have a heavy thermal protection systems. As shown here, 
an optimal design for the total vehicle weight may be obtained as discussed earlier. 
However, if the aeroassisted transfer is to be prefered to all propulsive transfer, it must 
offer a reduction in fuel mass greater than the increase in thermal protection mass. 
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TABLE 1 


EFFECTS OF HEAT CONSTRAINTS ON 
CHARACTERISTIC VELOCITIES 


PEAK HEATING RATES 

AV, 

AV 2 

AV 3 

TOTAL 

(WATTS/cm ) 

(m/sec) 

(m/sec) 

(m/s*c) 

(m/sec) 

UNCONSTRAINED 
6 = 239 

1031 

821 

18 

1870 

DR. NAIDU (BNDSCO) 
UNCONSTRAINED 

1034 

816 

43 

1893 

HEAT CONSTRAINT 
Q a 203 

1028 

855 

18 

1901 

HEAT CONSTRAINT 
6 = 170 

1026 

930 

19 

1974 



Fig. 1 Aeroassisted Orbital 
Plane Change 


TABLE 2 


CHARACTERISTIC VELOCITIES FOR 
DIFFERENT ORBIT INCLINATION CHANGES 


A| 

AVi 

AV 2 

ay, 

TOTAL AV 

(degree) 

(m/ sec) 

(m/sec) 

(m/sec) 

(m/sec) 

15 

1029 

312 

21 

1362 

20 

1031 

821 

18 

1871 

25 

1035 

1270 

31 

2336 



Fig. 2 Collocation and Hermite 
Approximation 
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